On the degeneracy of atomic states within exact-exchange (spin-) density functional theory 
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The problem of degenerate ground states of open-shell atoms is investigated in spin-restricted and 
unrestricted density functional theory using the exact exchange energy functional. For the spin- 
unrestricted case, spurious energy splittings of the order of 2 to 3 kcal/mol are found for atoms of 
the second and third period which is larger than the splittings obtained from recently proposed ap- 
proximate exchange functionals depending explicitly on the current density. In remarkable contrast, 
for spin-restricted calculations the degeneracy of different atomic ground states is recovered to within 
less than 0.6 kcal/mol. 



I. INTRODUCTION 

The Hohenberg-Kohn |1] and Kohn-Sham |2] theo- 
rems of density functional theory (DFT), which were 
originally established for non-degenerate ground states, 
may be extended to degenerate ground states as well 
00]. These degenerate ground states lead to a set of dif- 
ferent ground state densities and the exact energy func- 
tional yields the same ground state energy for all these 
densities. It has long been known, however, that com- 
mon approximations do not yield the same total ener- 
gies (al^, @]. In a systematic investigation of this 
problem Baerends and coworkers |9] showed that for 
states with different total magnetic quantum number, 
Ml, spurious energy splittings of up to 5 kcal/mol re- 
sult from generalized gradient approximations (GGA's). 
Even larger ones are observed for the meta-GGA's [10]. 

The problem has attracted renewed interest recently. 
Becke has proposed an approach for constructing ex- 
change-correlation functionals with an increased abil- 
ity to reproduce the degeneracy of atomic states filll . 
The essential idea is to enforce the proper description 
of the Fermi (or exchange) hole curvature 1 12] in the 
approximation of the exchange-correlation energy func- 
tional 1 13]. As a consequence, the paramagnetic current 
density appears explicitly in the expression of the corre- 
sponding functional 11411 . This improves the description 
of the atomic degeneracy ITlll of states carrying different 
paramagnetic current densities. 

Along the line of Becke's approach, Maximoff et al 
fl5ll have modified the GGA of Perdew, Burke, and 
Ernzerhof (PBE) to a form explicitly dependening on the 
current density. In this way they successfully reduced 
the previous spurious energy splittings. Actually, they 
have weakened Becke's suggestion by improving not 
the exchange hole curvature at all points in space, but 
rather its average. 

More recently, Tao and Perdew IToll have employed 
ideas of the current-DFT framework of Vignale and Ra- 
solt 1 16]. They constructed a current-dependent correc- 
tion to GGA and meta-GGA functionals and their results 
again suggest that some improvements for the energy 
splittings can be achieved. 

In this work we test the performance of the exact ex- 
change energy functional using the Optimized Effective 



Potential method |17, 18] for the problem of degener- 
ate ground states. In particular, we analyze an interest- 
ing aspect of the degeneracy problem related to the ad- 
ditional degrees of freedom introduced by going from 
the original DFT formulation of Hohenberg, Kobn and 
Sham to the spin-DFT (SDFT) formalism of von Barth 
and Hedin 11911 . As a consequence of this additional 
variational freedom, lower total energies are obtained in 
SDFT than in corresponding spin-restricted (DFT) calcu- 
lations. At the same time, however, the spurious energy 
splittings are increased for the states of different Ml- 
In this work we only consider densities represented by 
single Slater determinants of Kohn-Sham orbitals. The 
general formalism to deal with densities which can only 
be described as weighted sum of several determinantal 
densities is discussed in Refs. I20l l2lll . 

In the following we recall the basic ideas of the Opti- 
mized Effective Potential method and briefly compare 
the resulting equations in the DFT and SDFT frame- 
work. We then give some details on our numerical im- 
plementation along with the resulting energy splittings 
for the exact-exchange functional. Our findings are com- 
pared with results from other approximations discussed 
in the literature before we draw our conclusions. 



II. SPIN RESTRICTED AND UNRESTRICTED 
KOHN-SHAM SCHEMES 

In this Section we briefly review the basic equations of 
spin-restricted and spin-unrestricted density functional 
theory. We then focus on orbital-dependent approxi- 
mations to the exchange-correlation energy functional 
and discuss, for a given orbital functional, the relation 
between the corresponding exchange-correlation poten- 
tials in the restricted (DFT) and the unrestricted (SDFT) 
formalisms. 

In the unrestricted SDFT formalism, the total energy 
E of a system of interacting electrons is a functional of 
the two spin densities /v(r) (a =|, j): 



E iPhPl] = T s[phPl] 



d s r v (r)p(r) + U[p}+E xc [p hPl } 



(1) 



where 



T s [p h Pi\ =EE f d3r <PiA*) (-?) ^-w ( 2 ) 

is the non-interacting kinetic energy and N a is the num- 
ber of electrons with spin a. vq(t) is an external, electro- 
static potential and 



p{r) = p T (r) +p[(r) 



(3) 



is the total electronic density. The classical electrostatic 
(Hartree) interaction energy is given by 



U[p] 



d A r d 



„, P(r)p( r ') 



r - r' 



(4) 



and E xc is the exchange-correlation energy functional 
which has to be approximated in practice. The single- 
particle orbitals y>j„(r) in Eq. Q are solutions of the 
Kohn-Sham equation Ill9ll 



(5) 



where j is a collective index for the one-electron quan- 
tum numbers except spin. The effective single particle 
potential for spin a is given by 



v sa (r) = v (r) + v H (r) + v xca (v) 



with the Hartree potential 



v H (r) = jd 3 r 



, P(r') 



r — r 



and the exchange-correlation potential 

5E xc [p hPl ] 



5 pa{r) 



(6) 



(7) 



(8) 



The self-consistency cycle is closed by computing the 
spin densities via 



P<r( T ) = IViaW 

3=1 



(9) 



where the sum runs over the occupied orbitals. 

The unrestricted Kohn-Sham scheme of SDFT imme- 
diately reduces to the restricted scheme of DFT if one 
considers exchange-correlation functionals which only 
depend on the total electronic density of Eq. 10 as en- 
visioned in the original Hohenberg-Kohn theorem 1 1 ] . 
Then the exchange-correlation potential is 



5E xc [p] 
6p(r) 



(10) 



and both v xc and the total effective potential v s are in- 
dependent of the spin index a. Note that in spinte of 
the spin-independence of v s and v xc , the Kohn-Sham 
orbitals, being proper fermionic single-particle orbitals, 
still carry a spin-dependence. 

Of course, the exchange-correlation energy functional 
needs to be approximated in practice. Popular approxi- 
mations like the local density approximation (LDA) or 
generalized gradient approximations (GGA's) use an 
approximate form of E xc which explicitly depends on 
the density (in DFT) or on the spin-densities (in SDFT), 
i.e., different forms of the functional are required in DFT 
and SDFT. However, if one considers functionals which 
explicitly depend on the single-particle orbitals rather 
than the (spin-)density, one and the same orbital func- 
tional may be used either in the DFT or in the SDFT 
framework. The difference is the implicit dependence 
of the Kohn-Sham orbitals on the corresponding basic 
variables: in DFT they are implicit functionals of the to- 
tal particle density only, while in SDFT the orbitals are 
implicit functionals of the spin densities. 

For orbital functionals, the calculation of the ex- 
change-correlation potential is somewhat more com- 
plicated than for explicit density functionals and is 
achieved with the so-called Optimized Effective Poten- 
tial Method (OEP) E^IHl. For a review of the method 
the reader is referred to Ref. |22]. The OEP method leads 
to an integral equation for the exchange-correlation po- 
tential. For simplicity, we consider approximations of 
E xc that are functionals of the occupied orbitals only. 
The OEP integral equation can then be written in com- 
pact notation (in SDFT) 



£(^ CT (r)^(r)+ C .c.) =0. (11) 

3=1 

Here we have defined the orbital shifts 



^ CT (r) = /d 3 rV* CT (r') ( Vxca (r') - u xcj „(r')) G Sj(7 (r', r) 

(12) 

where Gsja is the Green function of the Kohn-Sham sys- 
tem 



G Sja (v',v)= J2 



k=l 



iL( r ')yv( r ) 



(13) 



and 



1 5E XC 



(14) 



In a series of steps 1 22, 23], the OEP equation can be 
transformed to 
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E 



\tp jcr {r)\ 2 (u xcja (r) + (vxcja- - u xcja )) - V • {ip* a (r)Vp ja (r)) 



+ c.c. 



(15) 



where 



and 



'Jxcja = Jd 3 r tp* a (r)v xccr (r)(p ja (r) , 



(16) 



xcja = \d 6 r ^ CT (r)u xcjo .(r)yi j(T (r) . 



(17) 



Similar expressions can, of course, be obtained for the 
spin-restricted case. The OEP equation analogous to 
Eq. ITU reads 



£ £(# ff (r)Mr)+cc.) =0 



(18) 



where the modified orbital shifts tpja are defined in anal- 
ogy to Eq. iT2i with being replaced by v xc . Apply- 
ing the same steps as in the SDFT case, the OEP equation 
of DFT transforms to 



v xc (r) 



2p(r) 



EE 

<x=U 3 = 1 



IWO)! 2 (u xc ja(r) + {Vxcja ~ Uxcja)) 



V (r)V<^(r)) 



+ c.c. 



(19) 



where v XC jc is defined as v xc j<r in Eq. fl6l except that 
v xca is again replaced by v xc . The DFT exchange- 
correlation potential <19t can be written as a weighted 
average of potentials for the different spin channels 



Pl(r)v xc i(r) + p 1 (r)f} a;ci (r) 
Pt( r ) + Pi( r ) 



(20) 



where 



2pa(r) 



E 



|w(r)| 2 (u XCJa (r) + (Vxcjv - Uxcje)) 



V-(^ CT (r)V^ ff (r)) 



(21) 



Eq. l2"0l shows how, in the spin-restricted case, the 
spin-up and spin-down channels mix to form the spin- 
independent exchange-correlation potential. 



Eqs. 11H or 1191 can be solved iteratively along 
with the corresponding Kohn-Sham equations in a self- 
consistent fashion. Due to the presence of the unoccu- 
pied Kohn-Sham orbitals in the definition of the orbital 
shifts (see Eqs. fi"2"l and fT3l), the full numerical solution 
of the OEP integral equation is nontrivial. In the origi- 
nal paper 1 18], solutions were presented for atomic sys- 
tems with spherical symmetry. Much later, it has also 
been solved for systems with lower symmetry such as 
molecules ] 24, 25] and solids |26]. Recently, an itera- 
tive algorithm for the solution of the OEP integral equa- 
tion based on the orbital shifts has been implemented 

In what follows we do not attempt a solution of the 
full OEP equation but rather use an approximation sug- 
gested by Krieger, Li and Iafrate |29] which has been 
found to be rather accurate in many situations. In this 
so-called KLI approximation, the terms containing the 
orbital shifts on the r.h.s. of Eqs. 1111 or <!19i are ne- 
glected completely. The KLI approximation may be sub- 
stituted by a slightly more elaborate one known as Com- 
mon Energy Denominator Approximation (CEDA) 1 3 1 
or Localized Hartree-Fock (LHF) approximation |3JJ 
However, it has been found that CEDA and KLI total 
energies are extremely close for atoms ] 32 ] . Moreover, 
for the atoms studied in this work we expect that KLI 
and CEDA results are very similar also for the current- 
carrying states since in most cases (from boron to mag- 
nesium) the current-carrying orbitals enter the expres- 
sions for the KLI and CEDA potentials in exactly the 
same way. 

The OEP equations given above are valid for any form 
of the exchange-correlation functional E xc which de- 
pends on the occupied orbitals only. In this work we 
use the exact exchange functional 



E x = - 



\H E 



<r=1,l3,k=l 



d A r d 



,3 ,^( r VL( r >^(r'W(r) 



(22) 



which is nothing but the Fock term of Hartree-Fock the- 
ory evaluated with Kohn-Sham orbitals. In our calcula- 
tions the correlation energy is neglected completely. 



III. NUMERICAL RESULTS 

In this Section we present our numerical results on the 
degeneracy problem of open-shell atomic ground states. 
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It is well-known that standard approximations like LDA 
or GGA do not give the same, degenerate ground state 
energies for different open-shell configurations. This is 
due to the fact that the densities of these ground states 
are different, leading to different Kohn-Sham potentials 
derived from these densities and therefore also to dif- 
ferent total energies [9]. Recently, this problem has at- 
tracted renewed interest |10, 11, 15] where approximate 
functionals depending on the current-density have been 
suggested which, while not solving the problem com- 
pletely, at least reduce the energy splittings between dif- 
ferent configurations significantly. Here we investigate 
the problem at the exact-exchange level, both in DFT 
and SDFT. 

Although the (interacting) Hamiltonian of an atom 
has spherical symmetry, the ground state densities of 
open-shell atoms typically are not spherical. However, 
for any of the possible degenerate ground states one can 
always find an axis for which the corresponding density 
exhibits cylindrical symmetry and we choose this axis 
as the z-axis of our coordinate system. We seek a Kohn- 
Sham single-particle potential with the same cylindrical 
symmetry. Then the magnetic quantum number to is a 
good quantum number to characterize the Kohn-Sham 
orbitals. We can then perform self -consistent calcula- 
tions by specifying how many orbitals with m = 0, 1, . . . 
be occupied for each spin channel and then keep this 
configuration fixed throughout the self-consistency cy- 
cle. For example, for the boron atom, one configuration 
has all spin-up electrons and the two spin-down elec- 
trons in to = states while in another configuration one 
of the spin-up electrons is required to occupy an rn = 1 
state with the other occupations unchanged. In this way 
current-carrying and zero-current states can be consid- 
ered. 

We have developed an atomic code for DFT and 
SDFT calculations in a basis set representation, assum- 
ing cylindrical symmetry of the Kohn-Sham potential. 
As basis functions we use Slater-type basis functions for 
the radial part multiplied with spherical harmonics for 
the angular part. We use the quadruple zeta basis sets 
(QZ4P) of Ref. |H for the Slater functions. 

We have tested our code by computing the total en- 
ergies of spherically symmetric atoms of the first and 
second row of the periodic table in exchange-only KLI 
approximation and compared with results from accu- 
rate, fully numerical codes available in the literature 
I22L I29I |3. Our code reproduces these energies to 
within a maximum deviation of 0.3 Kcal/ mol and an av- 
erage deviation of 0.1 Kcal /mol for the first-row atoms 
and to within a maximum deviation of 0.9 Kcal/mol and 
an average deviation of 0.5 Kcal/ mol for the second row. 
As a more relevant estimate of the accuracy of our cal- 
culations we have also computed the energy splittings 
between different configurations in LSD. Our results re- 
produce those reported in Ref. |11] with a deviation of 
less than 0.02 kcal/ mol. 

We then calculated self-consistent total energies for 



Atom A jB R " 


A b 
L\jPBE 


A SDFT 
^x-KLI 


A DFT 
^x-KLI 


B 0.62 


-0.4 


1.66 


0.06 


C 0.69 


-0.7 


1.58 


0.06 


O 1.23 


-0.6 


2.36 


0.55 


F 1.51 


-0.6 


2.32 


0.40 


Al 0.96 


0.2 


1.68 


0.04 


Si 0.84 


-0.1 


1.76 


0.05 


S 1.95 


0.7 


3.04 


0.34 


CI 1.73 


0.3 


3.15 


0.25 



"current-dependent exchange functional of Ref. flTH 
''current-dependent exchange functional of Eq.(17) of Ref. Il5ll 

TABLE I: Spurios energy splittings, A = E(\Ml\ = 1) — 
E(Ml = 0) in kcal/mol. All the columns, but the last one, 
refer to spin-unrestricted calculations 



the different configurations of open-shell atoms. Table 
|I] shows the energy differences (spurious energy split- 
tings) between Kohn-Sham Slater determinants with to- 
tal magnetic quantum number \M L \ = 1 and M L = 
in kcal/mol of our exchange-only KLI calculations of 
DFT and SDFT. For comparison we also list the re- 
sults of the current-dependent exchange-only function- 
als of Refs. (denoted jBR) and fil (denoted jPBE) 
in the first and second column, respectively. As can be 
seen, our SDFT results for the exact-exchange functional 
lead to larger splittings than both the jBR and the jPBE 
functionals. The idea behind the construction of these 
functionals is to improve the exchange-hole curvature 
by inclusion of the orbital paramagnetic current den- 
sity. Since in our calculations we have used the exact 
exchange functional (and therefore also the correct ex- 
change hole curvature) the success of the jBR and jPBE 
functionals in reducing the energy splittings might actu- 
ally be due to an overcorrection of their parent function- 
als. 

The most remarkable results of our calculations are 
the energy splittings obtained with a pure DFT (i.e., 
spin-restricted) calculation using the exact exchange 
functional (last column of Table0}. These spurious split- 
tings are in most cases more than an order of magnitude 
smaller than the corresponding SDFT results, therefore 
basically reproducing the exact degeneracy of the differ- 
ent ground-state configurations. Of course, due to the 
additional variational degree of freedom, total energies 
in SDFT are always lower than corresponding DFT re- 
sults. The price to be paid for this improvement, how- 
ever, are the unphysically increased energy splittings. 



IV. CONCLUSIONS 

In this work we have calculated the spurious energy 
splittings between atomic states of different quantum 
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number Ml using the exact-exchange functional. We 
have employed the KLI approximation to compute the 
Kohn-Sham exchange potential and found that for spin- 
unrestricted calculations the splittings are between 1 
and 3 kcal/mol for the atoms of the second and third 
period. Somewhat surprisingly, these are larger than 
the splittings reported with approximate exchange func- 
tionals which include the paramagnetic current density 
as an input parameter Illlll5ll . 

However, if the exact exchange functional is used in a 
spin-restricted DFT calculation, the spurious energy split- 
tings are reduced significantly, the largest one being of 
the order of 0.5 kcal/mol. One might speculate that the 



remaining splittings are due to the KLI approximation 
and could be further reduced if the full OEP equations 
for the exchange potential are solved. 
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